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Abstract. We present and compare three generically applicable signal processing 
. methods for periodic orbit quantization via harmonic inversion of semiclassical 

I recurrence functions. In a first step of each method, a band-limited decimated periodic 

04 ' orbit signal is obtained by analytical frequency windowing of the periodic orbit sum. In 

a second step, the frequencies and amplitudes of the decimated signal are determined 
by either Decimated Linear Predictor, Decimated Fade Approximant, or Decimated 
Signal Diagonalization. These techniques, which would have been numerically unstable 
^1^, without the windowing, provide numerically more accurate semiclassical spectra than 

' does the filter-diagonalization method. 

o 

^ ; PACS numbers: 05.45. -a, 03.65.Sq 

1. Introduction 

The semiclassical quantization of systems with an underlying chaotic classical dynamics 
is a nontrivial problem for the reason that Gutzwiller's trace formula |jl|, |^ does not 
usually converge in those regions where the eigenenergies or resonances are located. 
Various techniques have been developed to circumvent the convergence problem of 
periodic orbit theory. Examples are the cycle expansion technique [Q], the Riemann- 
Siegel type formula and pseudo-orbit expansions [Q, surface of section techniques [|], and 
a quantization rule based on a semiclassical approximation to the spectral staircase 0. 
These techniques have proven to be very efficient for systems with special properties, e.g., 
the cycle expansion for hyperbolic systems with an existing symbolic dynamics, while 
the other mentioned methods have been used for the calculation of bound spectra. 

Recently, an alternative method based upon Filter-Diagonalization (FD) has been 
introduced for the analytic continuation of the semiclassical trace formula [0, ^. The 
FD method requires knowledge of the periodic orbits up to a given maximum period 
(classical action), which depends on the mean density of states. The same holds true for 
the three methods presented in this paper. The semiclassical eigenenergies or resonances 
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are obtained by harmonic inversion of the periodic orbit recurrence signal. The FD 
method can be generally applied to both open and bound systems and has also proven 
powerful, e.g., for the calculation of semiclassical transition matrix elements and the 
quantization of systems with mixed regular-chaotic phase space [|1^]. For a review on 
periodic orbit quantization by harmonic inversion see |]Tl] . 

In this paper the techniques for the harmonic inversion of periodic orbit signals 
are further developed. The semiclassical signal, in action or time, corresponds to a 
"spectrum" or response in the frequency domain that is composed of a huge, in principle 
infinite, number of frequencies. To extract these frequencies and their corresponding 
amplitudes is a nontrivial task. In the previous work 0, |^, |ll| the periodic orbit signal 
has been harmonically inverted by means of FD [l^, |14| , which is designed for the 
analysis of time signals given on an equidistant grid. The periodic orbit recurrence signal 
is given as a sum over usually unevenly spaced 6 functions. A smooth signal, from which 
evenly spaced values can be read off, is obtained by a convolution of this sum with, e.g., 
a narrow Gaussian function. The disadvantages of this approach are twofold. Firstly, 
FD acts on this signal more or less like a "black box" and, as such, does not lend itself 
to a detailed understanding of semiclassical periodic orbit quantization. Secondly, the 
smoothed semiclassical signal usually consists of a huge number of data points. The 
handling of such large data sets, together with the smoothing, may lead to significant 
numerical errors in some of the results for the semiclassical eigenenergies and resonances. 

Here, we propose three alternative methods for the harmonic inversion of the 
periodic orbit recurrence signal that avoid these problems. In a first step we create 
a shortened signal which is constructed from the original signal and designed to be 
correct only in a window, i.e., a short frequency range of the total band width. Because 
the original signal is given as a periodic orbit sum of 6 functions, this "filtering" can be 
performed analytically resulting in a decimated periodic orbit signal with a relatively 
small number of equidistant grid points. In a second step the frequencies and amplitudes 
of the decimated signal are determined from a set of nonlinear equations. To solve 
the nonlinear system, we introduce three different processing methods: Decimated 
Linear Predictor (DLP), Decimated Fade Approximant (DPA), and Decimated Signal 
Diagonalization (DSD). The standard and well-known Linear Predictor (LP) and Fade 
Approximant (PA) would not have yielded numerically stable solutions if the signal 
had not first been decimated by the windowing (filtering) procedure. Furthermore, 
this separation of the harmonic inversion procedure into various steps may lead to a 
clearer picture of the periodic orbit quantization method. Numerical examples will 
demonstrate that the techniques proposed in this paper also provide more accurate 
results than previous applications of FD. 

The paper is organized as follows. In section ^ we briefly review the general idea of 
periodic orbit quantization by harmonic inversion. In section ^ we construct the band- 
limited decimated periodic orbit signal which is analyzed in section ^ with the help of 
three different methods, viz. DLP, DPA and DSD. In section |^ we present and compare 
results for the three-disk scattering system as a physical example and the zeros of the 
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Riemann zeta function as a mathematical model for periodic orbit quantization. Some 
concluding remarks are given in section ^. 

2. Periodic orbit quantization by harmonic inversion 

In order to understand the following, a brief recapitulation of the basic ideas of periodic 
orbit quantization by harmonic inversion is necessary. For further details see [|lT| . 

Following Gutzwiller fl], the semiclassical response function for chaotic systems 
is given by 

f%E)=g^,^{E) + J2A^oe''-° , (1) 

po 

where go'^{E) is a smooth function and S^o and are the classical actions and 
weights (including phase information given by the Maslov index) of the periodic orbit 
contributions. Equation (|l]) is also valid for integrable systems when the periodic orbit 
quantities are calculated not with Gutzwiller's trace formula, but with the Berry- Tabor 
formula [0 for periodic orbits on rational tori. The eigenenergies and resonances are 
the poles of the response function. Unfortunately, the semiclassical approximation (|1|) 
does not converge in the region of the poles and hence the problem is the analytic 
continuation of g'^'^{E) to this region. 

As done previously ^, we will also make the (weak) assumption that the 
classical system has a scaling property, i.e., the shape of periodic orbits does not depend 
on the scaling parameter, w, and the classical action scales as 

'S'po = wspo ■ (2) 

In scaling systems, the fluctuating part of the semiclassical response function, 

g^%w) = J2Apoe'-'-\ (3) 

po 

can be Fourier transformed readily to obtain the semiclassical trace of the propagator 



hoc 



C^'is) = ^ I g'\w)e-''^dw = A^oS {s - Spo) . (4) 



°° po 



The signal C^^{s) has 5-spikes at the positions of the classical periods (scaled actions) 
s = Spo of periodic orbits and with peak heights (recurrence strengths) ^po, i.e., C^'^{s) is 
Gutzwiller's periodic orbit recurrence function. Consider now the quantum mechanical 
counterparts of g^'^{w) and C^'^{s) taken as the sums over the poles Wk of the Green's 
function, 

g'^'^iw) = E — > (5) 

Cq-(s) = — / g'''^{w)e-'''"dw = -z V 4e-*""=' , (6) 

ZTT J-oo ^ 

with dk being the residues associated with the eigenvalues. In the case under study, 
i.e., density of state spectra, the d^ are the multiplicities of eigenvalues and should be 
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equal to 1 for non-degenerate states. Semiclassical eigenenergies Wk and residues dk 
can now, in principle, be obtained by adjusting the semiclassical signal, eq. (^), to the 
functional form of the quantum signal, eq. (j^), with the {wk,dk} being free generally 
complex frequencies and amplitudes. This procedure is known as "harmonic inversion" . 
The numerical procedure of harmonic inversion is a nontrivial task, especially if the 
number of frequencies in the signal is large {e.g., more than a thousand) or even infinite 
as is usually the case for periodic orbit quantization. Note that the conventional way 
to perform the spectral analysis, i.e., the Fourier transform of eq. (|) will bring us back 
to analyzing the non-convergent response function g^'^iw) in eq. (|^). The periodic orbit 
signal @ can be harmonically inverted by application of FD [[121 , |T3t |14| , which allows 
one to calculate a finite and relatively small set of frequencies and amplitudes in a given 
frequency window. The usual implementation of FD requires knowledge of the signal 
on an equidistant grid. The signal is not a continuous function. However, a smooth 
signal can be obtained by a convolution of C^^{s) with, e.g., a Gaussian function, 

C-{s) = -L- Y: Aoe"(^-^-)'/^'^' . (7) 
y/ Zira po 

As can easily be seen, the convolution results in a damping of the amplitudes, d^ — >■ 
di"^ = dk exp (-w;^cr^/2). The width a of the Gaussian function should be chosen 
sufficiently small to avoid an overly strong damping of amplitudes. To properly sample 
each Gaussian a dense grid with steps As ~ cr/3 is required. Therefore, the signal 
(0) analyzed by FD usually consists of a large number of data points. The numerical 
treatment of this large data set may suffer from rounding errors and loss of accuracy. 
Additionally, the "black box" type procedure of harmonic inversion by FD, which 
intertwines windowing and processing, does not provide any opportunity to seek a deeper 
understanding of semiclassical periodic orbit quantization. It is therefore desirable to 
separate the harmonic inversion procedure into two sequential steps: Firstly, the filtering 
procedure that does not require smoothing and, secondly, a procedure for extracting the 
frequencies and amplitudes. In section we will construct, by analytic filtering, a band- 
limited signal which consists of a relatively small number of frequencies. In section ^ 
we will present methods to extract the frequencies and amplitudes of such band-limited 
decimated signals. 



3. Construction of band-limited decimated signals by analytical filtering 



Consider a signal of a presumably large length N. We split the corresponding Fourier 
spectrum, which is also of length A^, into M frequency intervals. In general, a frequency 
filter can be applied to a given signal by application of the Fourier transform |16, p^ . 
Specifically, the signal, in time or action, is first transformed to the frequency domain, 
e.g., by application of the fast Fourier transform (FFT) or by using the closed-form 
expression, if available, for the Fourier integral. The transformed signal, which is 
essentially a low-resolution spectrum, is multiplied with a frequency filter function f{w) 
localized around a central frequency, wq, and zeroed out everywhere else outside the 
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selected window. This leads to a band-limited Fourier spectrum. The frequency filter 
f{w) can be rather general; typical examples are a rectangular window or a Gaussian 
function. The resulting filtered or band-limited spectrum is then shifted by —wq, 
relocated symmetrically around the frequency origin, w = 0, and transformed back 
to the time domain by the application of the inverse FFT to produce a band-limited 
signal valid only in the window defined by the filter function. Such a band-limited 
signal still has the same original length A^. It is at this step that we apply decimation 
which amounts to an enhancement of the sampling time by a factor of M thus leading 
finally to a band-limited decimated signal of length [N/M], where [u] denotes the integer 
part of u. In other words, since the band-width of the band-limited decimated signal 
is M times smaller than that of the original signal, the sampling rate, or dwell time, 
between signal samplings can now be M times larger. Hence, the filtered or band-limited 
signal can be reduced by retaining only those signal points for which the time or action 
indices are, say, 1, M + 1, 2M + 1, . . .. This technique is known as "beamspacing" 



or "decimation" |]T^ of band-limited signals. The flexibility in the choice of window size 
can ensure a numerically stable implementation of the processing methods presented 
below. 

The special form of the periodic orbit signal as a sum of S functions allows for 
an even simpler procedure, viz. analytical filtering. In the following we will apply a 
rectangular filter, i.e., f{w) = 1 for frequencies w G [wq — Aw, wq + Aw], and f{w) = 
outside the window. Generalization to other types of frequency filters is straightforward. 
Starting from the semiclassical response function (spectrum) g^'^{w) in eq. (^, which is 
itself a Fourier transform of the "signal" (^, and using a rectangular window we obtain, 
after evaluating the "second" Fourier transform, the band-limited (bl) periodic orbit 
signal, 

Cl^is) = — / g^'iw)e-''^---^^dw 



1 ^ rwo+Aw . 

27r pQ JuiQ— A«) 



^ ^^^^^ sin [{s-Spo)Aw] ^,,^^^^ ^ 

The introduction of Wq into the arguments of the exponential functions in (P), causes 
a shift of frequencies by —wq in the frequency domain. Note that Cbi(s) is a smooth 
function and can be easily evaluated on an arbitrary grid of points s„ < Smax provided 
the periodic orbit data are known for the set of orbits with classical action Spo < Smax- 
Applying now the same filter used for the semiclassical periodic orbit signal to the 
quantum one, we obtain the band-limited quantum signal 

1 nwo+Aw . 

C^'^is) = — / g^'^{w)e~''^'"""°^dw 

K 

= -lY. 4e-*(""=-"'»)^ , \wk-w^\<Aw. (9) 

k=l 
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In contrast to the signal C"^™(s) in eq. (|^), the band-hmited quantum signal consists of 
a finite number of frequencies Wk, k = 1, . . . , K, where K can be of the order of ~ (50- 
200) for an appropriately chosen frequency window, Aw. The problem of adjusting 
the band-limited semiclassical signal in eq. (j^) to its quantum mechanical analogue in 
eq. can now be written as a set of 2K nonlinear equations 

C-(nr) = On = -1^2 4e-™>^ , n = 0,1, . . . ,2K - 1 , (10) 

k=l 

for the 2K unknown variables, viz. the shifted frequencies, w'l^ = Wk—Wo, and amplitudes, 
dk- The signal now becomes "short" (decimated) as it can be evaluated on an equidistant 
grid, s = UT, with step width r = tt/Aw. It is important to note that the number of 
signal points c„ in eq. (p!OD is usually much smaller than a reasonable discretization 
of the signal Cl^{s) in eq. (|^, which is the starting point for harmonic inversion by 
FD. Therefore, the discrete signal points c„ = C^i(nr) are called the "band-limited 
decimated" periodic orbit signal. Methods to solve the nonhnear system, eq. (|T(]|), are 
discussed in section || below. 

It should also be noted that the analytical filtering in eq. (§) is not restricted to 
periodic orbit signals, but can be applied, in general, to any signal given as a sum of 



S functions. An example is the high resolution analysis of quantum spectra [^, [11 
where the density of states reads g{E) = J2n ^{E — E. 



-'nj 



4. Harmonic inversion of decimated signals 

In this section we want to solve the nonlinear set of equations 

K 

On = Y.dkzl, n = 0,l,...,2K-l, (11) 

k=l 

where Zk = exp (— iw^r) and dk are generally complex variational parameters. For 
notational simplicity we have absorbed the factor of — z on the right-hand side of eq. (p!0|) 
into the dkS with the understanding that this should be corrected for at the end of the 
calculation. We assume that the number of frequencies in the signal is relatively small 
{K ~ 50 to 200). Although the system of nonlinear equations is, in general, still 
ill-conditioned, frequency filtering reduces the number of signal points, and hence the 
number of equations. Several numerical techniques, that otherwise would be numerically 
unstable, can now be applied successfully. In the following we introduce three different 
methods, viz. Decimated Linear Predictor (DLP), Decimated Fade Approximant (DFA), 
and Decimated Signal Diagonalization (DSD). 



4.1. Decimated Linear Predictor 

The problem of solving eq. ( pA]) has already been addressed in the 18th century by 
Baron de Frony, who converted the nonlinear set of equations (|TI|) to a linear algebra 
problem. Today this method is known as Linear Predictor (LP). Our method, called 
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Decimated Linear Predictor (DLP), strictly applies the procedure of LP except with one 
essential difference; the original signal C^'^{s) is replaced by its band-limited decimated 
counterpart c„ = (nr). 

Equation ( [Tl| ) can be written in matrix form for the signal points Cn+i to Cn+K, 



( c„+i \ 



V Cn+K j 



i 



^1 



\ ^1 



( \ 

V dK ] 



(12) 



From the matrix representation (|T2D it follows that 



/ 



\ ^1 



-K 





( Cn+1 \ 


) 


\ Cn+K / 



K 

CtkCn+k 

k=l 



(13) 



which means that every signal point c„ can be "predicted" by a linear combination of 
the K subsequent points with a fixed set of coefficients Ok, k = 1, . . . , K. The first step 
of the LP method is to calculate these coefficients. Writing eq. (0) in matrix form with 
n = 0,...,K — l,we obtain the coefficients Ok as solution of the linear set of equations, 



\ Ck 



ck \ 



/ ai \ 



( '° ^ 

\ Ck-1 J 



(14) 



■ C2K-1 J \ CtK J 

The second step is the determination of the parameters Zk in eq. (pT]). Using eqs. (0) 
and ( |Tl| ) we obtain 

K K K 

Cn = J2 (^kCn+k = ^Y1 (^kdlZ^j"^^ , (15) 



and thus 



K 

E 

fc=i 



k=l 



K 



1=1 k=l 



.1=1 



.,n+l 



dk = 0. 



(16) 



Equation (|T6D is satisfied for arbitrary sets of amplitudes dk when Zk is a zero of the 
polynomial 



K 

1=1 



. 



:i7) 



The parameters Zk = exp (— zw^r) and thus the frequencies 



log(2fc) 



are therefore obtained by searching for the zeros of the polynomial in eq. ([171) . Note 
that this is the only nonlinear step of the algorithm and numerical routines for finding 
the roots of polynomials are well established. In the third and final step, the amplitudes 
dk are obtained from the linear set of equations 



K 



E 44 



n 



0, 



(19) 



k=l 
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To summarize, the LP method reduces the nonlinear set of equations (0) for the 
variational parameters {zk,dk} to two well-known problems, i.e., the solution of two 
linear sets of equations (0) and ([T9| ) and the root search of a polynomial, eq. ([T7|), 
which is a nonlinear but familiar problem. The matrices in eqs. ( ]T^ ) and ( |I9| ) are 
a Toeplitz and Vandermonde matrix, respectively, and special algorithms are known 
for the fast solution of such linear systems [T^\. However, when the matrices are ill- 
conditioned, conventional LU decomposition of the matrices is numerically more stable, 
and, furthermore, an iterative improvement of the solution can significantly reduce errors 
arising from numerical round-off. The roots of polynomials can be found, in principle, 
by application of Laguerre's method [l^. However, it turns out that an alternative 
method, i.e., the diagonalization of the Hessenberg matrix 

/ 



V 



a-K-i 
aK 
1 









a.K-2 
aK 









ai 




aK 


aK 














1 


J 


) = 


detiA 



(20) 



for which the characteristic polynomial P{z) = det[A — zl] =0 is given by eq. (p!7|) 
(with Oo = —1), is a numerically more robust technique for finding the roots of high 
degree {K > 60) polynomials |jl9[ . 



4.2. Decimated Pade Approximant 

As an alternative method for solving the nonlinear system ([TI|) we now propose to 
apply the method of Decimated Pade Approximants (DPA). This is the standard Pade 
Approximant (PA) but applied to our band-limited decimated signal c„. Let us assume 
for the moment that the signal points c„ are known up to infinity, n = 0,1,... 00. 
Interpreting the the coefficients of a Maclaurin series in the variable 2; ^, we 

can then define the function g{z] 
geometric series we obtain 



J2n=0 ^nZ 



K 



With eq. ( |Tl|) and the sum rule for 

_ Pk{z) 



K 



zdk 



n=0 fc=l n=0 k=l ^ ' 



Qk{z) 



(21) 



The right-hand side of eq. (^) is a rational function with polynomials of degree K in 
the numerator and denominator. Evidently, the parameters Zk = exp (— zty^r) are the 
poles of g{z), i.e., the zeros of the polynomial Qk{z). The parameters dk are calculated 
via the residues of the last two terms of (pi]) . We obtain 

ZkQ'xiz^k) 

with the prime indicating the derivative d/dz. Of course, the assumption that the 
coefficients c„ are known up to infinity is not fulfilled and, therefore, the sum on the 
left-hand side of eq. (^l]) cannot be evaluated in practice. However, the convergence of 
the sum can be accelerated by application of DPA. Indeed, with DPA, knowledge of 2K 
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signal points Cq, . . . , C2K-1 is sufficient for the calculation of the coefficients of the two 
polynomials 

K K 

Pk{z) = E hz"" and Qk{z) = ^ " 1 • (23) 

fc=l k=l 

The coefficients Ok-, k = 1, . . . ,K are obtained as solutions of the linear set of equations 

K 

Cn = E "-kCn+k , U = 0, . . . , K - 1 , 

k=l 

which is identical to eqs. (0) and (|14|) for DLP. Once the a's are known, the coefficients 
bk are given by the explicit formula 

K-k 

bk= (^k+mCm , k = 1, . . . ,K . (24) 

m=0 

It should be noted that the different derivations of DLP and DPA provide the same 
polynomial whose zeros are the Zk parameters, i.e., the calculated with both methods 
exactly agree. However, DLP and DPA do differ in the way the amplitudes, dk, are 
calculated. It is also important to note that DPA is applied here as a method for signal 
processing, i.e., in a different context to that in ref. pOf, where the Pade approximant is 



used for the direct summation of the periodic orbit terms in Gutzwiller's trace formula. 
4.3. Decimated Signal Diagonalization 

In refs. []T2|, |I^] it has been shown how the problem of solving the nonlinear set of 
equations (|TT|) can be recast in the form of the generalized eigenvalue problem, 

\JBu = zuSBu . (25) 

The elements of the K x K operator matrix U and overlap matrix S depend trivially 
upon the c„'s |T^ : 



Uij — Cj+j+i ; Sij — Ci^j ; i, j — 0, . . . , K — 1 . (26) 

Note that the operator matrix U is the same as in the linear system ([T^), i.e. the matrix 
form of eq. (|13D of DLP. The matrices U and S in eq. (^) are complex symmetric [i.e., 
non-Hermitian), and the eigenvectors are orthogonal with respect to the overlap 
matrix S, 

iBk\S\Bk,) = NAk' , (27) 

where the brackets define a complex symmetric inner product {a\b) = {b\a), i.e., no 
complex conjugation of either a or b. The overlap matrix S is not usually positive 
definite and therefore the N^s are, in general, complex normalization parameters. An 
eigenvector Bk cannot be normalized for Nk = 0. The amplitudes dk in eq. (|11]) are 
obtained from the eigenvectors Bk via 

\-K-l 

(28) 



dk 

Nk 



1 

^ ^ ^n-Bk,n 
n=0 
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The parameters Zk in eq. (|TT]) are given as the eigenvalues of the generahzed eigenvalue 
problem (^), and are simply related to the frequencies w'j. in eq. (|1^) via = 
exp(-iw^r). 

The three methods introduced above (DLP, DPA and DSD) look technically quite 
different. With DLP the coefficients of the characteristic polynomial ( [T7| ) and the 
amplitudes dk are obtained by solving two linear sets of equations (0) and (|TP|). Note 
that the complete set of zeros z^. of eq. ([T7| ) is required to solve for the d^ in eq. ([19|) . 
The DPA method is even simpler, as only one linear system, eq. (0), has to be solved 
to determine the coefficients of the rational function Pk{z)/Qk{z). Finding the zeros 
of eq. ( |17|) gives knowledge about selected parameters Zk, and allows one to calculate 
the corresponding amplitudes dk via eq. (^). The DSD method requires the most 
numerical effort, because the solution of the generalized eigenvalue problem (|25|) for 
both the eigenvalues Zk and eigenvectors Bk is needed. 

It is important to note that the three methods, in spite of their different derivations, 
are mathematically equivalent and provide the same results for the parameters {zk, dk}, 
when the following two conditions are fulfilled: the nonlinear set of equations (|Tl|) has 
a unique solution, when, firstly, the matrices U and S in eq. ( p6D have a non-vanishing 
determinant (det U 7^ 0, det S 7^ 0), and, secondly, the parameters Zk are non-degenerate 
{zk 7^ Zk' for k 7^ k'). These conditions guarantee the existence of a unique solution 
of the linear equations (p!^ ) and (p!9|), the non-singularity of the generalized eigenvalue 
problem (pSf), and the non- vanishing of both the derivatives Qxi^k) in eq. (p2D and the 
normalization constants Nk in eqs. (pT]) and (pSl). Equation (|1T|) usually has no solution 
in the case of degenerate Zk parameters. However, degeneracies can be handled with 
a generalization of the ansatz (|TlD and modified equations for the calculation of the 
parameters. This special case will be reported elsewhere. 

While the parameters Zk in eq. (|Tll) are usually unique, the calculation of the 
frequencies w'^. via eq. (|18D is not unique, because of the multivalued property of the 
complex logarithm. To obtain the "correct" frequencies it is necessary to appropriately 
adjust the range Aw of the frequency filter and the step width r of the band-limited 
decimated signal (plQl). We recommend the following procedure. The most convenient 
approach is to choose first the centre Wq of the frequency window and the number K 
of frequencies within that window. Note that K determines the dimension of the linear 
systems, and hence the degree of the polynomials which have to be handled numerically, 
and is therefore directly related to the computational effort required. Frequency windows 
are selected to be sufficiently narrow to yield values for the rank between if ^ 50 and 
K ^ 200. The step width for the decimated signal should be chosen as 



with Sjnax being the total length of the periodic orbit signal. The relation Zk = 
exp {—iw'i^T) projects the frequency window w' G [—Aw, +Aw] onto the unit circle 
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in the complex plane when the range of the frequency window is chosen as 

TT 27rK , , 

Aw = - = . (30) 

When calculating the complex logarithm with arglogz G [— vr, +7r], eq. ([TSD provides 
the "correct" shifted frequencies w'^. and thus the frequencies Wk = Wi^\ w'^. 

To achieve convergence, the length Smax of the periodic orbit signal must be 
sufficiently long to ensure that the number of semiclassical eigenvalues within the 
frequency window is less than K. As a consequence the harmonic inversion procedure 
usually provides not only the true semiclassical eigenvalues but also some spurious 
resonances. The spurious resonances are identified by low or near zero values of the 
corresponding amplitudes d^ and can also be detected by analyzing the shifted decimated 
signal, i.e., signal points ci, . . . , c^k instead of cq, . . . , C2k-\- The true frequencies usually 
agree to very high precision, while spurious frequencies show by orders of magnitude 
larger deviations. 

5. Results and discussion 

In this section we want to demonstrate the efficiency and accuracy of the method 
introduced above by way of two examples: the three-disk repeller as an open physical 
system and the zeros of the Riemann zeta function as a mathematical model for 
periodic orbit quantization of bound chaotic systems. Both systems have previously 
been investigated by means of FD 0, |^, [ll|], which allows us to make a direct comparison 



of the results. The three-disk scattering system has also served as a prototype for the 
development and application of cycle expansion techniques [0, ^ , and we will briefly 



discuss the differences between harmonic inversion and cycle expansion. 
5.1. The three-disk repeller 

As the first example, we consider a billiard system consisting of three identical hard 
disks with unit radii, = 1, displaced from each other by the same distance d. This 
simple, albeit nontrivial, scattering system has already been used as a model within 
the cycle expansion method |P, 21, ^ and periodic orbit quantization by harmonic 



inversion [|^, |TT|. We give therefore only a very brief introduction to the system 
and refer the reader to the literature for details. The three-disk scattering system is 
invariant under the symmetry operations of the group Cs^,, i.e., three refiections at 
symmetry lines and two rotations by 27r/3 and 47r/3. Resonances belong to one of the 
three irreducible subspaces A\, A2, and E p3[. As in most previous work we concentrate 
on the resonances of the subspace Ai for the three-disk repeller with R = 1 and d = 6. 

In billiards, which are scaling systems, the shape of the periodic orbits does not 
depend on the energy E, and the classical action is given by the length L of the orbit, i.e., 
Spo = wspo = hkLpo (see eq. (^), where w = k = \k\ = ylME/% is the absolute value 
of the wave vector to be quantized. Setting ^ = 1, we use Spo = Lpo in what follows. 
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In figure [l|a we present the periodic orbit signal C^'^{L) for tlie tliree-disk repeller in tlie 
region < L < Lmax = 35. Tlie signal is given as a periodic orbit sum of delta functions 
5{L — Lpo) weighted with the periodic orbit amplitudes (see eq. (||)). The groups 
with oscillating sign belong to periodic orbits with adjacent cycle lengths. Signals of this 
type have been analyzed (after convolution with a narrow Gaussian function, see eq. (|^) 
by FD in refs. 0, ^, ^, 0, |TT|. We now illustrate harmonic inversion of band- limited 
decimated periodic orbit signals by DLP, DPA and DSD. 

In a first step, a band-limited decimated periodic orbit signal is constructed as 
described in section ^ As an example we choose K = 100 as the rank of the nonlinear 
set of equations (|T0|) , and ko = 200 as the centre of the frequency window. The width of 
the frequency window is given by Ak = 2'nK/L^a_^ = 2007r/35 ~ 18.0. The step width 
of the decimated signal is r = AL = Lmax/2-ft' = 0.175. The band-limited decimated 
periodic orbit signal points Cn = Cli{L = nAL), with n = 0, . . . , 2K are calculated with 
the help of eq. (|]) and presented in figure |l|b. The solid and dashed lines are the real 
and imaginary parts of C^f(?7.AL), respectively. The modulations with spacings n/Ak 
result from the superposition of the sinc-like functions in eq. (P). 

The band-limited decimated periodic orbit signal C^i(nAL) can now be analyzed, 
in a second step, with one of the harmonic inversion techniques introduced in section ^, 
viz. DLP, DPA or DSD. The resonances obtained by DLP are presented as plus symbols 
in figure |T]c. The dotted lines at Re = 182 and Re = 218 show the borders of 
the frequency window. The two symbols very close to the border on the left indicate 
spurious resonances. 

A long range spectrum can be obtained by choosing several values Wq for the 
centre of the frequency window in such a way that the windows slightly overlap. 
Figure ^ presents the semiclassical resonances for the three-disk repeller in the range 
< Re < 250. The spectrum has been obtained by harmonic inversion of decimated 
periodic orbit signals similar to that in figure |l|b but with an increased signal length, 
Lmax = 55. The plus symbols, crosses, and squares denote the semiclassical resonances 
obtained by DLP, DPA and DSD, respectively. The resonances obtained by the three 
different harmonic inversion techniques are in perfect agreement, with the exception of 
a few resonances in the region Re /c < 25. In this region the matrices U and S in 
eq. ( p6D are rather ill-conditioned, and the few discrepancies can therefore be explained 
as numerical artifacts. 

The spectrum presented in figure was obtained previously in ref. [|| by application 
of FD |12, 13, 14|. In table |l] we compare the semiclassical eigenvalues k and residues 
dk of selected resonances obtained by (a) FD and (b) harmonic inversion of band- 
limited decimated periodic orbit signals. For non-degenerate resonances under study 
the residues should he d^ = I. In [|] the residues of several resonances deviate 
significantly from dk = 1 hj more than 5% (see the resonances marked by (") in 
table 0). With harmonic inversion of decimated signals the accuracy of the residues is 
increased by several orders of magnitude (see the resonances marked by (^) in table |I[). 
The semiclassical eigenvalues k also reveal deviations between the different numerical 
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techniques. The resonances obtained by harmonic inversion of band-hmited decimated 
signals are in much better agreement with the results obtained by the cycle expansion 



method [22| than those obtained by FD in ref. 

Numerical values for the residues very close to dk = 1 indicate well converged 
semiclassical resonances, and this is the case for all resonances of the four bands closest 
to the real axis in figure |^. Resonances with non-vanishing dk have also been obtained 
in the region Re > 120, Im k < —0.8 (see figure H). These resonances, although not 



fully converged, are in qualitative agreement with exact quantum calculations p2[. It 
is important to note that the different techniques for harmonic inversion of decimated 
signals, viz. DLP, DPA and DSD, yield the same results, even for those resonances 
which are not fully converged. This illustrates the mathematical equivalence of the 
three methods as explained in section |. 

5.2. Harmonic inversion vs. cycle expansion 

The three-disk scattering system discussed in section 0| has purely hyperbolic classical 



dynamics and has been used extensively as the prototype model within the cycle 
expansion techniques |P, As has been shown by Voros Gutzwiller's trace 



formula for unstable periodic orbits can be recast in the form of an infinite and non- 
convergent Euler product over all periodic orbits. When the periodic orbits obey a 
symbolic dynamics the semiclassical eigenenergies or resonances can be obtained as 
the zeros of the cycle expanded Gutzwiller- Voros zeta function. Unfortunately, the 
convergence of the cycle expansion is restricted, due to poles of the Gutzwiller- Voros zeta 



function . The domain of analyticity of semiclassical zeta functions can be extended 



25| , p6| resulting in the "quasiclassical zeta function" [^, ^], which is an entire function 



for the three-disk repeller. This approach allows one to calculate semiclassical resonances 
in critical regions where the Gutzwiller- Voros zeta function does not converge, at the 
cost, however, of many extra spurious resonances and with the rate of convergence 
slowed down tremendously p2 . 



With the limited numerical accuracy of harmonic inversion by FD applied in ref. |^ , 
the semiclassical resonances of the three-disk repeller in the region Im < —0.6 were 
somewhat unreliable. The improved accuracy of the analysis of band- limited decimated 
periodic orbit signals introduced in the present paper now allows us to compare the 
two semiclassical quantization techniques, viz. harmonic inversion and cycle expansion 
methods, even for resonances deep in the complex plane. We will demonstrate that the 
harmonic inversion method provides semiclassical resonances in energy regions where 
the cycle expansion of the Gutzwiller- Voros zeta function does not converge. 

In figure § we present a part of the semiclassical resonance spectrum of figure § in 
the region 25 < Re A; < 65. The squares and crosses label the semiclassical resonances 
obtained by harmonic inversion of the decimated semiclassical periodic orbit signal and 
cycle expansion of the Gutzwiller- Voros zeta function |2^ , respectively. The dotted line 
in figure ^ indicates the borderline, Im k = —0.121 557 0, which separates the domain 
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of absolute convergence of Gutzwiller's trace formula from the region where analytic 
continuation is necessary. For the two resonance bands slightly below this border the 
results of both semiclassical quantization methods are in perfect agreement. The dashed 
line in figure marks the abscissa of absolute convergence for the Gutzwiller-Voros 



zeta function at Im /c = —0.699110 |2^. The Gutzwiller-Voros zeta function provides 
several spurious resonances which accumulate at Im /c ~ —0.9, i.e., slightly below the 
borderline of absolute convergence (see the crosses in figure H). The resonances in the 
region Im k < —0.9, especially those belonging to the fourth band, are not described by 
the Gutzwiller-Voros zeta function but are obtained by the harmonic inversion method 
(see the squares in figure 1^). 



5.3. Zeros of the Riemann zeta function 

As the second example to demonstrate the numerical accuracy of harmonic inversion of 
band-limited decimated periodic orbit signals we investigate the Riemann zeta function 
which is a mathematical model for periodic orbit quantization. Here we only briefly 
explain the idea of this model and refer the reader to the literature 



rT| for details. 



The hypothesis of Riemann is that all the nontrivial zeros of the analytic 
continuation of the function 



c(^) = E 



n 



n(i 



p 



n=l 

have real part ^, so that the values w 



-1 



2' 



(Re z > 1, p : primes) 
Wk, defined by 



c 







(31) 



(32) 



are all real or purely imaginary ^8 
poles of the function 



where 



P m=l 



\og{p) 



The parameters Wk can be obtained as the 



(33) 



p 



m/2 



^pm 



(34) 
(35) 



mlog(p) , 

with p indicating the prime numbers. As was already pointed out by Berry |2£ 
eq. ( P^ has the same mathematical form as Gutzwiller's trace formula with the primes 
interpreted as the primitive periodic orbits, Apm and Spm the "amplitudes" and "classical 
actions" of the periodic orbit contributions, and m formally counting the "repetitions" of 
orbits. Equation (^) converges only for Imw > ^ and analytic continuation is necessary 
to extract the poles of g{w), i.e., the Riemann zeros. The advantage of studying the 
zeta function instead of a "real" physical bound system is that no extensive periodic 
orbit search is necessary for the calculation of Riemann zeros, as the only input data are 
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just prime numbers. Harmonic inversion can be applied to adjust the Fourier transform 
of eq. ([33|), i.e., 

oo 

Cis) =J2Y1 ^pmSis - Spm) , (36) 
p m=l 

to the functional form 

where the are the Riemann zeros and the residues d^ have been introduced as 
adjustable parameters which here should all be equal to 1. 

In ref. about 2600 Riemann zeros to about 12 digit precision have been obtained 
by harmonic inversion of the signal (^) with Smax = log(lO^) ~ 13.82 using FD. 
However, the numerical residues agree with dk = 1 only to about a 5 or 6 digit precision. 
With harmonic inversion of band-limited decimated signals the accuracy of both the 
Riemann zeros Wk and the multiplicities dk is improved by several orders of magnitude. 
In table | we compare selected values obtained by (a) FD |§] and (b) DLP. The increase 
in accuracy can easily be seen, in particular for the imaginary parts, Im Wk and Im dk, 
which should both be equal to zero. The same improvement in accuracy is also achieved 
by application of DPA and DSD. 

The precise calculation of parameters dk = I for the residues of the Riemann zeros 
does not seem to be of great interest. However, it should be noted that the multiplicities 
may be greater than one, e.g., for some eigenvalues of integrable systems such as the 
circle billiard ||lT], ^ where states with angular momentum quantum number m ^ are 
twofold degenerate. As we have checked the techniques presented in this paper indeed 
yield the correct multiplicities to very high precision. The dkS have also nontrivial 
values when used, e.g., for the semiclassical calculation of diagonal matrix elements 
and non-diagonal transition strengths M in dynamical systems. 



6. Conclusion 



We have introduced three methods for semiclassical periodic orbit quantization, 
viz. Decimated Linear Predictor (DLP), Decimated Pade Approximant (DPA), and 
Decimated Signal Diagonalization (DSD) for the harmonic inversion of band-limited 
decimated periodic orbit signals. The characteristic feature of these methods is the strict 
separation of the two steps, viz. the analytical filtering of the periodic orbit signal and 
the numerical harmonic inversion of the band-limited decimated signal. The separation 
of these two steps and the handling of small amounts of data compared to other "black 
box" type signal processing techniques enables an easier and deeper understanding of 
the semiclassical quantization method. Furthermore, applications to the three-disk 
repeller and the Riemann zeta function demonstrate that the new methods provide 
numerically more accurate results than previous applications of filter-diagonalization 
(FD). A detailed comparison of various semiclassical quantization methods reveals that 
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quantization by harmonic inversion of the band-hmited decimated periodic orbit signal 
can even be apphed in energy regions where the cycle expansion of the Gutzwiller-Voros 
zeta function does not converge. 

The methods introduced in this paper can be applied to the periodic orbit 
quantization of systems with both chaotic and regular classical dynamics, when the 
periodic orbit signal is calculated with Gutzwiller's trace formula |^ for isolated orbits 
and the Berry- Tabor formula JTSf for orbits on rational tori, respectively. More generally, 
any signal given as a sum of S functions can be filtered analytically and analyzed using 
the methods described in sections ^ and ^. For example, the technique can also be 
applied to the harmonic inversion of the density of states g{E) = J2n^{E — En) of 
quantum spectra to extract information about the underlying classical dynamics [|l^, [Tl[]. 

It is to be noted that all the signal processing techniques mentioned in this paper 
lend themselves to a formulation where non-diagonal responses appear in the frequency 
domain and their complementary cross-correlation-type signals appear in the time or 
action domain ||3l|, |3^. This is important as such methods allow for the use of shorter 
signals and hence, in the context of this paper, fewer periodic orbits. The need to find 
large numbers of periodic orbits may limit the practical utility of these methods and so 
any attempt to overcome this problem is worth investigating. Cross-correlation methods 
also sample better and yield improved results for poles (resonances) that lie deep in the 
complex plane than do straight correlation-based signal processing techniques. 
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Figure 1. (a) Periodic orbit recurrence signal for the three-disk scattering system 
with R = 1, d = 6 without filtering. The signal in the region i < 35 consists of 93 

non-equidistant periodic orbit contributions (including multiple repetitions), (b) Same 
as (a) filtered with frequency window w G [182,218]. The decimated signal consists 
of 201 equidistant data points with AL = 0.175. The solid and dashed lines are the 
real and imaginary part of C{L), respectively, (c) Semiclassical resonances obtained 
by harmonic inversion of the decimated signal C{L) in (b). The dotted lines mark the 
borders of the frequency window. 
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Figure 2. Semiclassical resonances for the three-disk scattering system {Ai subspace) 
with R = 1, d = Q obtained by harmonic inversion via Decimated Linear Predictor 

(phis symbols), Decimated Pade Approximant (crosses), and Decimated Signal 
Diagonahzation (squares) of the analyticaUy decimated periodic orbit signal. 
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Figure 3. Semiclassical resonances {Ai subspace) for the three-disk scattering system 
witli R = I, d = 6. Squares: Harmonic inversion of the decimated semiclassical 
recurrence signal; Crosses: Cycle expansion of the Gutzwiller-Voros zeta function 
[ p2[ . The dotted and dashed lines mark the borderline for absolute convergence of 
Gutzwiller's trace formula (Im k — —0.121 557) and the Gutzwiller-Voros zeta function 
(Im k = —0.699110), respectively. The harmonic inversion method converges deeper 
in the complex plane than the Gutzwiller-Voros zeta function. 
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Table 1. Semiclassical resonances and multiplicities for the three-disk scattering 
problem (Ai subspace) with R ~ I, d — 6. (°) Filter-diagonalization method (FD) 
(^) Decimated Linear Predictor (DLP). 
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Table 2. Nontrivial zeros Wk and multiplicities for the Riemann zeta fmiction. 
(") Filter-diagonalization method (FD) [^; (^) Decimated Linear Predictor (DLP). 
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